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Abstract 

We discuss how to generate random unitary matrices from the classical com- 
pact groups U(A^), O(A^) and USp(2A^) with probability distributions given by the 
respective invariant measures. The algorithm is straightforward to implement us- 
ing standard linear algebra packages. This approach extends to the Dyson circular 
ensembles too. This article is based on a lecture given by the author at the sum- 
mer school on Number Theory and Random Matrix Theory held at the University 
of Rochester in June 2006. The exposition is addressed to a general mathematical 
audience. 



1 Introduction 

Since Random Matrix Theory (RMT) was introduced by Wishart [17J in 1928, it has found 
applications in a variety of areas of physics, pure and applied mathematics, probability, 
statistics and engineering. Few examples — far from being exhaustive — include: analytic 
number theory, combinatorics, graph theory, multivariate statistics, nuclear physics, quan- 
tum chaos, quantum information, statistical mechanics, structural dynamics and wireless 
telecommunications. The main reasons for the ever growing success of RMT are mainly 
two. Firstly, in the limit of large matrix dimension the statistical correlations of the spectra 
of a family, or ensemble, of matrices are independent of the probability distribution that 
defines the ensemble, but depend only on the invariant properties of such a distribution. 
As a consequence random matrices turn out to be very accurate models for a large num- 
ber of mathematical and physical problems. Secondly, RMT techniques allow analytical 
computations to an extent that is often impossible to achieve in the contexts that they are 
modelling. This predictive ability of RMT is particularly powerful whenever in the original 
problem there are not any natural parameters to average over. 

Although the advantage of using RMT lies in the possibility of computing explicit math- 
ematical and physical quantities analytically, it is sometimes necessary to resort to numer- 
ical simulations. The purpose of this article is twofold. Firstly, we provide the reader with 
a simple method for generating random matrices from the classical compact groups that 
most mathematicians — not necessarily familiar with computer programming — should 
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be able to implement in a code of only few lines. This is achieved in section [5l Secondly, 
we discuss in detail the main ideas, which turn out be fairly general and quite interesting, 
behind this algorithm. 

An N X N unitary matrix U = {ujk) is defined by the relation U*U = UU* = I, which 
in terms of the matrix elements reads 

N N N N 

^u*f,Uki = ^UkjUki = 5ji and ^Ujkuli = ^UjkUik = 5ji, (1.1) 

k=l k=l k=l k=l 

where U* = (Wj^) is the conjugate transpose of U, i.e. u*^ = Ukj- In this article we will 
use the symbol "to denote complex conjugation, in order to distinguish it from *, which is 
reserved to the conjugate transpose of a matrix. The constraints (11.11) simply state that 
the columns (rows) of a unitary matrix form an orthonormal basis in C^. The set U(A^) of 
unitary matrices forms a compact Lie group whose real dimension is A^^; it is then made 
into a probability space by assigning as a distribution the unique measure invariant under 
group multiplication, known as Haar measure. Such a probability space is often referred 
to as Circular Unitary Ensemble (CUE). 

Usually the correct ensemble to model a particular situation depends on the symmetries 
of the problem. Ensembles of unitary matrices are constructed in two steps: we first identify 
a subset U C U(A^) by imposing further restrictions on U; then we assign to U a probability 
measure with the appropriate invariant properties. As well as U(A^), we will discuss how 
to generate random matrices from the orthogonal O(A^) and unitary symplectic USp(2A^) 
groups with probability distributions given by the respective unique invariant measures. 
We shall also consider the two remaining Dyson circular ensembles [2], namely the Circular 
Orthogonal Ensemble (COE) and Circular Symplectic Ensemble (CSE). Other symmetric 
spaces appear in the applications [18], but we will not concern ourselves with them. 

Writing an algorithm to generate random unitary matrices which is both correct and 
numerically stable presents some pitfalls. The reason is that the conditions (11. ip imply 
that the matrix elements are not independent and thus are statistically correlated. The 
main ideas discussed in this article are centred around the QR decomposition and go back 
to Wedderburn [16], Heiberger [5] (corrected by Tanner and Thisted [IS]), Stewart [S] 
and Diaconis and Shahshahani fL\. However, the technical literature may be difficult to 
access for a reader without a background in numerical analysis or statistics, while the im- 
plementation of such techniques is elementary. Another method discussed in the literature 
involves an explicit representation of the matrix elements of U in terms of A^^ independent 
parameters (Euler angles) [19], but it does not seem to be equally efficient or convenient. 

2 Some examples and motivations 

Before discussing how to generate random matrices it is helpful to give few examples that 
show how they appear in the applications. 
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In quantum mechanics all the information about an isolated physical system at a given 
time to is contained in a state vector ipo belonging to a Hilbert space H — in general 
infinite dimensional. The time evolution of tpo, i-^- its dynamics, is determined by a 
unitary operator U. In other words, at a time t > to, i/jq has evolved into 

^ = U^o. (2.1) 

The fact that U is unitary guarantees that W^pW = Wi'oW = 1, which is required by the 
probabilistic interpretation of quantum mechanics. 

If the dynamics is complicated — as in heavy nuclei or in quantum systems whose clas- 
sical limits are characterized by a chaotic dynamics — writing down an explicit expression 
for U may be hopeless. Therefore, we can attempt to replace [/ by a random operator and 
check if the predictions that we obtain are consistent with the empirical observations. It 
is also reasonable to simplify the problem even further and replace [/ by a random unitary 
matrix of finite, but large, dimension. Then the main question is: What are the matrix 
space and the probability distribution that best model our system? 

In physics the symmetries of a problem are often known a priori, even if the details 
of the dynamics remain obscure. Now, suppose that our system is invariant under time 
reversal but does not have any other symmetry. From general considerations (see Mehta [H] 
p. 36) we know that U is always conjugate by a unitary transformation to a symmetric 
matrix. Therefore, we can always choose U so that 

U = U\ (2.2) 

where [/* denotes the transpose of U. Since there are not other symmetries in the problem, 
this is the only constraint that we can impose. Therefore, the appropriate matrices that 
model this physical system should be symmetric. Let us denote by the set of unitary 
symmeytric matrices. If t/ e it can be proved (see Meta [9] p. 499) that it admits the 
representation 

U = WW\ W e UiN). (2.3) 

This factorization is not unique. Let O(A^) be the group of real matrices O such that 
00* = 0*0 = / and set W = WO. By definition we have 

U = W'W'^ = WOO^W* = WW\ (2.4) 

This statement is true also in the opposite direction: if WW* = WW* there exists an 
O G O(A^) such that W' = WO. Therefore, is isomorphic to the left coset space of O(A^) 
in U(iV), i.e. 

= U(A^)/0(A^). (2.5) 

Since a measure space with total mass equal to one is a probability space, in what follows 
we shall use the two terminologies interchangeably. An ensemble of random matrices is 
defined by a matrix space and a probability measure on it. We have found the former; we 
are left to identify the latter. Haar measure, which will be discussed in detail in section [SI 
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provides a natural probability distribution on U(A^); 'natural' in the sense that it equally 
weighs different regions of U(A^), thus it behaves like a uniform distribution. From the 
factorization (12. 3p the probability distribution on U(A^) induces a measure on 0. As a 
consequence, if W is Haar distributed the resulting measure on will be uniform too. In 
section [S] we shall see that such a measure is the unique probability distribution induced 
by Haar measure on 0. Therefore, it provides a natural choice to model a time reversal 
invariant quantum system. The space together with this measure is the COE ensemble. 

If a quantum system does not have any symmetry, then there are no restriction to 
U(A^) and the natural choice of probability distribution is Haar measure. This is the CUE 
ensemble. If the system is invariant under time reversal and has a half-integer spin, then 
the appropriate ensemble is the CSE. The matrix space of the CSE is the subset S C U(2A^) 
whose elements admit the representation 

U = -WJW'J, W G U(2A^), (2.6) 

where 

From the factorization (12. 6p the probability distribution on U(2iV) induces a measure on 
S. As previously, such a measure is fixed by assigning Haar measure to U(2iV). 

The set S is isomorphic to a coset space too. The unitary symplectic group USp(2A^) 
is the subgroup of U(2A^) whose elements obey the relation 

SJS' = J. (2.8) 

Therefore, the matrix U in equation (12. 6p does not change if we replace W with W' = WS, 
where S G USp(2A^). Similarly, if W and W are such that 

U = -WJW'J = -W'JW'J, W, W G U(2A^), (2.9) 

then W'W~^ G USp(2iV). Therefore, 

S ^ U(2A^)/USp(2Ar). (2.10) 

The probability distribution of the CSE is the unique invariant measure induced on the 
coset space (I2.10p by Haar measure on U(2A^). 

From equations (12. 3p and (12. 6p all we need to generate random matrices in the CUE, 
COE and CSE ensembles is an algorithm whose output is Haar distributed unitary matrices. 
The rest of this article will concentrate on generating random matrices from all three 
classical compact groups U(A^), O(A^) and USp(2A^) with probability distributions given 
by the respective Haar measures. These groups are not only functional to constructing 
matrices in the COE and CSE, but are also important ensembles in their own right. Indeed, 
the work of Montgomery [TT], Odlyzko [12], Katz and Sarnak [6], Keating and Snaith [ZIIH] 
and Rubinstein has shown beyond doubt that the local statistical properties of the 
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the Riemann zeta function and other L-functions can be modelled by the characteristic 
polynomials of Haar distributed random matrices. Over the last few years the predictive 
power of this approach has brought about an impressive progress in analytic number theory 
that could not have been achieved with traditional techniques. (See [TU] for a collection of 
review articles in the subject.) 



3 Haar measure and invariance 



Since the algorithm the we shall discuss is essentially based on the invariant properties of 
Haar measure, in this section we introduce the main concepts that are needed to understand 
how it works. We, nevertheless, begin with another ensemble: the Ginihre ensemble. 
Besides being a simpler illustration of the ideas we need, generating a matrix in the Ginibre 
ensemble is the first step toward producing a random unitary matrix. 

The space of matrices for the Ginibre ensemble is GL(iV, C), the set of all the invert- 
ible N X N complex matrices Z = {zjk)\ the matrix elements are independent identically 
distributed (i.i.d.) standard normal complex random variables. In other words, the prob- 
ability density function (p.d.f.) of zjk is 

piz.k) = -e~l^-"=r (3.1) 

71 

By definition the matrix entries are statistically independent, therefore the joint probability 
density function {j.p.d.f.) for the matrix elements is 

1 ^ 2 1 / ^ \ 1 

^(^) = ^ n = ^^^P - E l^^-^l' = ^exp(-TrZ*Z). (3.2) 

j,k=l \ j,k=l / 

Since P{Z) is a probability density, it is normahzed to one, i.e. 

P{Z)dZ=l, (3.3) 

where dZ = Y[fk=i'^^jkdyjk and zjk = xjk + iyjk- The j.p.d.f. P{Z) contains all the 
statistical information on the Ginibre ensemble. 

Since C^^^ = , we will use the two notations according to what is more appropriate 
for the context. Thus, we can write 



d^iG{Z) = P{Z)dZ (3.4) 



and think of rf/iQ as an infinitesimal volume or measure in 
we say that d^Q is invariant under / if 



d^^G{f{z))=d^lG{z). 



(3.5) 



6 



Francesco Mezzadri 



Lemma 1. The measure of the Ginibre ensemble is invariant under left and right multi- 
plication of Z by arbitrary unitary matrices, i.e. 

dficiUZ) = dficiZV) = dficiZ), U,Ve U(iV). (3.6) 

Proof. First we need to show that P{UZ) = P{Z); then we must prove that the Jacobian 
of the map 

Z^UZ (3.7) 
(seen as a transformation in C^^) is one. Since by definition U*U = I, we have 

PiUZ) = exp (- Tr Z*U*UZ) = exp (- Tr Z*Z) = P{Z). (3.8) 

Now, the map fl3.7p is isomorphic to 

X = U ®---®U. (3.9) 

A''times 

It follows immediately that X is a N"^ x N'^ unitary matrix, therefore |detX| = 1. The 
proof of right invariance is identical. □ 

Because the elements of a unitary matrix are not independent, writing an explicit for- 
mula for the infinitesimal volume element of U(X) is more complicated than for the Ginibre 
ensemble. An N x N unitary matrix contains 2N'^ real numbers and the constraints (11. ip 
form a system of N'^ real equations. Therefore, U(X) is isomorphic to a iV^-dimensional 
manifold embedded in M^^^. Such a manifold is compact and has a natural group structure 
that comes from matrix multiplication. Thus, an infinitesimal volume element on U{N) 
will have the form 

djjiiU) = m(ai, . . . , aj\f2)dai ■ ■ ■ daN2, (3.10) 

where ai, . . . ,aN2 are local coordinates on U(X). Every compact Lie group has a unique 
(up to an arbitrary constant) left and right invariant measure, known as Haar measure. In 
other words, if we denote Haar measure on U(X) by dfiu^U), we have 

duniVU) = dfiuiUW) = dfiuiU), V,W E \]{N). (3.11) 

Although an explicit expression for Haar measure on U(X) in terms of local coordinates 
can be written down (see Zyczkowski and Kus [19j for a formula), we will see that in order 
to generate matrices distributed with Haar measure we only need to know that is invariant 
and unique. 

Haar measure normalized to one is a natural choice for a probability measure on a 
compact group because, being invariant under group multiplication, any region of U(X) 
carries the same weight in a group average. It is the analogue of the uniform density on 
a finite interval. In order to understand this point consider the simplest example: U(l). 
It is the set {e'^} of the complex numbers with modulo one, therefore it has the topology 
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of the unit circle Since in this case matrix multiphcation is simply addition mod 27r, 
U(l) is isomorphic to the group of translations on S^. A probability density function that 
equally weighs any part of the unit circle is the constant density p{9) = l/(27r). This is 
the standard Lebesgue measure, which is invariant under translations. Therefore, it is the 
unique Haar measure on U(l). 

Note that it is not possible to define an 'unbiased' measure on a non-compact manifold. 
For example, we can provide a finite interval with a constant p.d.f. p{x), but not the whole 
real line M, since the integral p{x)dx would diverge. 

4 The QR decomposition and a numerical experiment 

By definition the columns of a x unitary matrix are orthonormal vectors in C^. Thus, 
if we take an arbitrary complex N x N matrix Z of full rank and apply the Gram-Schmidt 
orthonormalization to its columns, the resulting matrix Q is unitary. It turns out that if 
the entries of Z ares i.i.d. standard complex normal random variables, i.e. if Z belongs 
to the Ginibre ensemble, then Q is distributed with Haar measure (see Eaton [S], p. 234, 
for a proof). Unfortunately, the implementation of this algorithm is numerically unstable. 
However, we may observe that 

Z = QR, (4.1) 

where R is upper-triangular and invertible. In other words, the Gram-Schmidt algorithm 
realizes the QR decomposition. This factorization is widely used in numerical analysis to 
solve linear least squares problems and as first step of a particular eigenvalue algorithm. 
Indeed, every linear algebra package has a routine that implements it. In most cases, 
however, the algorithm adopted is not the Gram-Schmidt orthonormalization but uses the 
Householder reflections, which are numerically stable. 

Because of this simple observation, at first one might be tempted to produce a matrix 
in the Ginibre ensemble and then to use a black box QR decomposition routine. Writing 
such a code is straightforward. For example, if we choose the SciPy library in Python, we 
may implement the following function: 

from scipy import * 

def wrong_distribution(n) : 

>'''>' Random matrix with the wrong distribution' ' ' ' ' ' 

z = (randn(n,n) + lj*randn(n,n))/sqrt(2.0) 

q,r = linalg.qr(z) 

return q 

Unfortunately, the output is not distributed with Haar measure, as it was observed by 
Edelman and Rao [1]. It is instructive to give an explicit example of this phenomenon. 

A unitary matrix can always be diagonalized in U(A^). Therefore, its eigenvalues 
{e*^^, . . . , e*^^} lie on the unit circle. A classical calculation in RMT (see Mehta P] pp. 203- 
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205) consists of computing the statistical correlations among the arguments 6j. The sim- 
plest correlation function to determine is the density of the eigenvalues p{0), or — as 
sometimes it is called — the one-point correlation. Since Haar measure is the analogue of 
a uniform distribution, each set of eigenvalues must have the same weight, therefore the 
normalized eigenvalue density is 

Pi0) = ^- (4.2) 

ZTT 

It is important to point out that equation (14.21) does not mean that the eigenvalues are 
statistically uncorrelated. 

Testing (14.21) numerically is very simple. We generated 10, 000 random unitary matrices 
using wrong_distribution(n). The density of the eigenvalues of such matrices is clearly 
not constant (figure 1(a)). Figure 1(b) shows the histogram of the spacing distribution, 
which deviates from the theoretical prediction too. This statistics is often plotted because it 
encodes the knowledge of all the spectral correlations and is easy to determine empirically. 
For unitary matrices it is defined as follows. Take the arguments of the eigenvalues and 
order them in ascending order: 

Oi<e2<...<eN. (4.3) 
The normalised distances, or spacings, between consecutive eigenvalues are 

N 

^. = ^(^.+1-^.), J = h...,N. (4.4) 

The spacing distribution p{s) is the probability density of s. (For a discussion on the 
spacing distribution see Mehta j3| p. 118.) 

It is worth emphasising the QR decomposition is a standard routine. The most com- 
monly known mathematical software packages like Matlab, Mathematica, Maple, SciPy for 
Python essentially use a combination of algorithms found in LAPACK routines. Changing 
software would not alter the outcome of this numerical experiment. 



5 A correct and efficient algorithm 

What is wrong with standard QR factorization routines? Where do they differ from the 
Gram-Schmidt orthonormalization? Why is the probability distribution of the output 
matrix not Haar measure? 

The main problem is that QR decomposition is not unique. Indeed, let Z G GL(A^, C) 
and suppose that Z = QR, where Q is unitary and R is invertible and upper-triangular. If 



A 



\ 



diag(e'^\...,e*^^) 



(5.1) 
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(a) Eigenvalue density (b) Spacing distribution 

Figure 1: Empirical histograms of the density of the eigenvalues and of the spacing distri- 
butions compared with the theoretical curves for the CUE. The data are computed from 
the eigenvalues of ten thousand 50 x 50 random unitary matrices obtained from the routine 
wrong_distribution(n) . 



then 

Q' = QA and R' = A'^R (5.2) 
are still unitary and upper-triangular respectively. Furthermore, 

Z = QR = Q'R'. (5.3) 

Therefore, the QR decomposition defines a multi-valued map 

QR:GL(iV,C) — >V{N)xT{N), (5.4) 

where T(A^) denotes the group of invertible upper-triangular matrices. 

In order to make the mapping (15.41) single-valued, we need to specify the algorithm 
that achieves the factorization. In most applications such a choice is dictated only by the 
performance and stability of the code. For our purposes, however, the subset of 
U(iV) X T(A^), in which the output of the QR decomposition is chosen, is fundamental 
and we need to pay particular attention to it. It is convenient from a mathematical point 
of view to introduce a variation of the mapping (15.41) . which is not only single- valued but 
also one-to-one. In this way we will not have to refer all the times to a specific algorithm. 
Indeed, the idea is that we should be able to alter the output of a QR decomposition 
routine without even knowing the algorithm implemented. 

We first need 

Lemma 2. Equation (15.31) implies (15. 2p . where A G A(A^) and A{N) denotes the group of 
all unitary diagonal matrices (15.11) . 
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Proof. Equation (15.31) can be rearranged as 

g-ig' = RR'-\ (5.5) 

Since U(A^) and T(A^) are groups, both sides of equations (15. 5p must belong to U(A^)nT(A^). 
By definition the inverse of a unitary matrix U is its conjugate transpose and the inverse 
of an upper-triangular matrix is upper-triangular. Therefore, if a matrix is both unitary 
and upper-triangular it must be diagonal, i.e. A{N) = U(A^) fl T(A^). □ 

This lemma suggests that, more naturally, instead of the QR factorization (15. 4p we 
should consider a one-to-one map 

QR:GL(A^,C) — >V{N)xT{N), (5.6) 

where r(A^) = T(A^)/A(A^) is the right coset space of A{N) in T{N). We construct (I5J]) as 
follows: we first define it on a class of representatives of T{N) using the QR factorization; 
then we extend it to the whole r(A^). However, since the QR decomposition is not unique, 
there is a certain degree of arbitrariness in this definition. We need to find a map under 
which the measure of the Ginibre ensemble induces Haar measure on U(A^). The main tool 
to achieve this goal is the invariance under group multiplication of Haar measure and its 
uniqueness. Thus, our choice of the decomposition (15.61) must be such that if 

Z^(g,7) then UZ ^ (UQ,^) (5.7) 

with the same 7 G r(A^) and for any U G U(A^). This property implies that left multi- 
plication of Z by a unitary matrix reduces, after the decomposition, to the left action of 
U(A^) into itself. But lemma [H states that 

d^^o{UZ) = dficiZ) (5.8) 

for any U G U(A^). As a consequence, if the map (15. 6p satisfies (15. 7p the induced measure on 
U(A^) will be invariant under left multiplication too and therefore must be Haar measure. 

How do we construct the map (15.61) ? A class of representatives of r(A^) can be chosen 
by fixing the arguments of the elements of the main diagonal of i? G T(A^). Let us impose 
that such elements all be real and strictly positive. Using (15. 2p we can uniquely factorize 
any Z G GL(A^, C) so that the main diagonal of R has this property. It follows that if 
Z = QR, then 

UZ = UQR, U G U(A^). (5.9) 

This QR decomposition of UZ is unique within the chosen class of representatives of r(A^). 
Therefore, the resulting map (15.61) obeys (15. 7p . Finally, we arrive at 

Theorem 1. Suppose that the map (15.61) satisfies the hypothesis (15.71) . Then, it decomposes 
the measure (13.40 of the Ginibre ensemble as 



duG^Z) = dfiuiQ) X dfir{N){l)- 



(5.10) 
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Proof. We have 

duoiU Z) = diiQ^Z) by lemma [H (5.11a) 

= dfi{UQ,'y) = dfi{Q,'y) by equation (15. 7p (5.11b) 

= dnn{Q) X dfir(N){l) by the uniqueness of Haar measure. (5.11c) 

□ 

The choice of the class of representatives that we made coincides exactly with outcome 
of the Gram-Schmidt orthonormalization. The output of standard QR decomposition 
routines are such that if Z ^ {Q, R) then UZ ^ {Q\ R') with Q' 7^ UQ and R' ^ R. 
Therefore, the corresponding map (15.61) does not obey (15. 7p and theorem [1] does not hold. 

We can now give a recipe to create a random unitary matrix with distribution given by 
Haar measure. 

1. Take an NxN complex matrix Z whose entries are complex standard normal random 
variables. 

2. Feed Z into any QR decomposition routine. Let {Q,R), where Z = QR, be the 
output. 

3. Create the following diagonal matrix 

I \rii\ \ 



(5.12) 



A = 

\ kivjvl/ 
where the T jj s circ the diagonal elements of R. 

4. The diagonal elements of R' = A~^R are always real and strictly positive, therefore 
the matrix Q' = QA is distributed with Haar measure. 

The corresponding Python function is: 

from scipy import * 
def haar_measure (n) : 

>'>'>> Random matrix distributed with Haar measure' ' ' > > ' 

z = (randn(n,n) + lj*randn(n,n))/sqrt(2.0) 

q,r = linalg.qr(z) 

d = diagonal(r) 

ph = d/absolute(d) 

q = multiply(q,ph,q) 

return q 

If we repeat the numerical experiment discussed in section H] using this routine, we 
obtain the histograms in figure 2, which are consistent with the theoretical predictions. 
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(a) Eigenvalue density (b) Spacing distribution 

Figure 2: Empirical histograms of the density of the eigenvalues and of the spacing distri- 
butions compared with the theoretical curves for the CUE. The data are computed from 
the eigenvalues of ten thousand 50 x 50 random unitary matrices output of the function 
haar_measure (n) . 

6 The unitary symplectic group USp(2A^) 

Up to now we have only considered U(A^). The discussion for O(A^) is identical, except that 
the input matrix of the QR decomposition routine must be real. Unfortunately, however, 
for USp(2A^) there are not any black box routines that we can use and we must put more 
effort into writing an algorithm. 

The algebra of unitary symplectic matrices can be rephrased in terms of Hamilton's 
quaternions; it is convenient for our purposes to use this formalism. A quaternion g G EI 
is a linear combination 

q = a ■ 1 + bii + ci2 + di^, a, b,c,d & M, (6.1) 
where 1 is the identity and zi, Z2, ^3 are the quaternion units; they obey the algebra 

= ^2 = "^3 = '^i^2^3 = —1- (6.2) 
We can also define the conjugate of q, 

q = a ■ 1 — bii — ci2 — di^, (6.3) 

as well the norm 

\\qf = qq = qq = a"^ + b'^ + + d'^, (6.4) 

When c = d = 0, M reduces to C and q is simply the complex conjugate of q. 

In analogy with and — provided we are careful with the fact that multiplication 
in EI is not commutative — we can study the space H^. Elements in EI^ are A^-tuples 
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q = (gi, . . . , qn)- The bilinear map 

N 

{p,^) = J2m^ p,qeH^, (6.5) 

is the analogue of the usual Hermitian inner product in and the norm of a quaternion 
vector is simply 

N 

l|qf = (q,q) = 5^ Ikif • (6.6) 

Similarly, GL(A^, H) is the group of all the N x N invertible matrices with quaternion 
elements. 

The quaternion units admit a representation in terms of the 2x2 matrices 

/2=(J °), ei=(j ''={-1 J) ^""^ '=^=0 o)' ^^•'^^^ 

where 

1 I— >• I2, ii I— >• ei, 12 I— >■ 62 and 1— >• 63. (6.7b) 
Thus, q = a ■ 1 + bii + ci2 + di^ is mapped into the complex matrix 




(6.8a) 



where z = a + ih and w = c + id. In addition 

Equations (16.81) generalize to an arbitrary N x N quaternion matrix Q, which can be 
represented in terms of a 2N x 2N complex matrix Q using the decomposition 

Q ^ Q = Qo ® ^2 + <5i ® ei + Q2 ® 62 + Qs ® eg, (6.9) 

where Qo, Qi, Q2 and Qs are arbitrary N x N real matrices. Proceeding in the same 
fashion, if Q G GL(A^, H) we define its conjugate transpose Q* = {q*f^) by setting g*^ = g^.^-. 

The symplectic group Sp(A^) is the subset of GL(iV, H) whose matrices satisfy the 
identity S*S = SS* = I. Because of the analogy between U(A^) and Sp(A^), the latter is 
sometimes called hyper-unitary group and is denoted by U(A^, H). The usefulness of the 
quaternion algebra lies in 

Theorem 2. The groups Sp(A^) andUSp{2N) are isomorphic, i.e. 

USp(2iV) = Sp(A^). (6.10) 
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Proof. It is convenient to replace the skew-symmetric matrix J in the definition (12.81) with 

/O 1 \ 

-1 ■•• 



n 



V 



1 
-1 0/ 



/® 62. 



(6.ii: 



This substitution is equivalent to a permutation of the rows and columns of J, therefore 
it is simply a conjugation by a unitary matrix. 

We first prove that if 5 G Sp(A^), then its complex representation S belongs to 
USp(2A^). By equation fl6.9l) S* is mapped to 



S^Q® h- Si ^ei- 3^2^62- Si® eg 
which follows from the identities 



-QS'Q, 



(A (g) 5)* = A* ® 5* and {A ® B){C D) = AC ® BD, 



(6.12) 



(6.13) 



and from the algebra (16.21) of the quaternion units. As a consequence, 

SS* ^ -SQS'Q = I. (6.14) 
Therefore, the matrix S is symplectic. Combining equations (]6.8bl) and (16.121) gives 

-QS'Q = S*, (6.15) 

thus S e USp(2A^). 

We now need to show that if S* G USp(2A^) then it is the representation of a matrix 
S G Sp(A^). This statement follows if we prove that S admits a decomposition of the 
form (16.91) . where 5*0, Si, S2 and S3 N x N must be real matrices. If this is true, then the 
argument of the first part of the proof can simply be reversed. 

Let us allow the coefficients a, b, c and d in the definition (16.11) to be complex numbers. 
The definitions of conjugate quaternion (16. 3p and conjugate transpose of a quaternion 
matrix, however, remain the same. The matrices (I6.7al) form a basis in C^^^. Therefore, 
any 2x2 complex matrix can be represented as a linear combination of I2, ei, 62 and 63. 
Thus, any matrix Q G C^^^^^ admits a decomposition of the form (16. 9p . but now the 
matrices Qo, Qi, Q2 and Q3 are allowed to be complex. In other words, Q is always the 
representation of a quaternion matrix Q, but in general the quaternion units have complex 
coefficients. The important fact that we need to pay attention to is that 



(6.16) 



if and only if the coefficients of the quaternion units are real numbers. This is a straight- 
forward consequence of the representation (16. Sal) . 
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Let S G USp(2A^) be the complex representation of the quaternion matrix S, but 
assume that S* is not mapped into S*. It is still true, however, that 

S*^-QS'Q, (6.17) 

because equation 06.121) is only a consequence of matrix manipulations. But since S is 
unitary symplectic S* = —Q Q, which is a contradiction. □ 

The algebra of Sp(A^) is the generalization to Hamilton's quaternions of the algebra 
of U(A^). Therefore, it is not surprising that the discussion of section [5] is not affected by 
replacing GL{N,C) and U{N) with GL(A^,H) and Sp(A^) respectively. Thus, since Sp(A^) 
and USp(2A^) are isomorphic, USp(2A^) and Sp(A^) have the same Haar measure dfiu- In 
particular, we can introduce the quaternion Ginibre ensemble, which is the set GL(A^, H) 
equipped with the probability density 

P{2) = -^exp(-TrZ*Z) = exp (- \\z,,\A . (6.18) 

\ j,k=l / 

Quaternion matrices can be factorized by the QR decomposition too: for any 
Z G GL(A^, H) we can always write 

Z = QTZ, (6.19) 

where Q G Sp(A^) and TZ is an invertible and upper-triangular quaternion matrix. Now, 
let 

A(iV,e) = {AGT(iV,e)|A = diag(gi,...,g^), ||g,|| = 1, j = l,...,Ar}, (6.20) 

where T(A^, H) is the group of invertible upper-triangular quaternion matrices. Further- 
more, let r(A^,H) = T(A^,H)/A(A^,M) be the right coset space of A(A^,e) in T(A^,H). We 
have the following 

Theorem 3. There exists a one-to-one map 

QTZ : GL(iV, e) — > Sp(iV) X r(iV, H) (6.21) 

such that 

Z^{Q,-f) and UZ ^ {UQ,^), (6.22) 

where 7 = K{N ,M)TZ. Furthermore, it factorizes the measure dfiQ of the Ginibre ensemble 
as 

dficiZ) = dfiiiiQ) X rf/ir(7v,e)(7)- (6-23) 
We leave proving these generalizations as an exercise for the reader. 
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7 Householder reflections 



Theorem [3] provides us with the theoretical tools to generate a random matrix in USp(2A^). 
However, when we implement these results in computer code, we need to devise an algo- 
rithm whose output satisfies the condition (16.221) . The first one that comes to one's mind is 
the Gram-Schmidt orthonormalization. But given that black box routines for quaternion 
matrices do not exists on the market, and that we are forced to write the complete code 
ourselves, we may as well choose one that is numerically stable and which, as it turns out, 
requires the same effort. The most common algorithm that achieves the QR decomposition 
uses the Householder reflections. For the sake of clarity, we will discuss this method for 
O(A^); the generalizations to U(A^) and Sp(A^) are straightforward. 

Given an arbitrary vector v G M"*, the main idea of the Householder reflections is to 
construct a simple orthogonal transformation Hm (dependent on v) such that 



|v|| ei. 



(7.1) 



where ei = (1,0,..., 0) G M™. For any real matrix X = (xjk), Hn is determined by 
replacing v in equation (17.11) with the flrst column of X. The product H^X will have the 
structure 

/rii * ■ ■ ■ *\ 
* ■ ■ • * 



HnX 



(7.2) 



where 



Then, deflne the matrix 



where 



Hn-1 




Hn-i{^'W 



V ei. 



(7.3) 



(7.4) 



(7.5) 



In this case v' is the (A^ — l)-dimensional vector obtained by dropping the flrst element of 
the second column of the matrix (17. 2p . We proceed in this fashion until the matrix 



R = H1H2 ■ ■ ■ Hn~iHnX 



(7.6) 



is upper-triangular with diagonal entries rn, . . . , tnn- The product 

/O Tjt fjt fjt fjt 



(7.7) 
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X2 




is by construction an orthogonal matrix. In equations 07.61) and (17.71) Hm denotes the block 
matrix 



m 



where Hm is defined in equation (17. ip . 

The matrix Hm is constructed using elementary geometry. Consider a vector in the 
plane, v = {xi,X2), and assume, for simplicity, that xi > 0. Furthermore, let u denote the 
unit vector along the interior bisector of the angle that v makes with the Xi-axis, i.e. 

^ v+ ||v|| ei 

u=M — — u, (7.9) 

||v + ||v|| ei|| 

where ei is the unit vector along the xi-axis. The refiection of v along the direction 
of u is ||v|| ei (see figure 3). Refiections are distance-preserving linear transformations, 
therefore their representations in an orthonormal basis are orthogonal matrices. In this 
simple example it can be constructed from elementary linear algebra: 

i/2(v) = -/ + 2uu*. (7.10) 

Finally, we obtain 

/72(v)v = ||v|| ei. (7.11) 

It is worth noting that H2 depends only on the direction of v and not on its modulus. 
Thus we can rewrite equation (17.111) as 



i72(v)v = ei. 



(7.12) 
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where v = v/ ||v||. 

The generahzation to arbitrary dimensions is straightforward. For any vector v G 
the Householder reflection is defined as 



where 



Furthermore, we have 



^n^(v) = T(/-2uU% (7.13) 



V ± ||v|| ei 

a= II , II - (7.14) 

V ± V ei 



Hm{^)v = ei. (7.15) 

How do we choose the sign in the right-hand side of equation fl7.14p ? From a mathe- 
matical point of view such a choice is irrelevant: in both cases Hmi'v) maps v into a vector 
whose only component different from zero is the first one. However, numerically it can be 
important. The square of the denominator in (17.141) is 

||v± ||v||ei||^ = 2||v|| (||v|| ±a;i), (7.16) 

where xi is the first component of v. If xi is comparable in magnitude to ||v|| and negative 
(positive) and we choose the plus (minus) sign, then the term 

||v||±a;i, (7.17) 

could very small and cancellations with significant round-off errors may occur. Therefore, 
the Householder transformation to be implemented in computer code should be 

H^iv) = - sgn(xi) (/ - 2uu*) , (7.18) 

where 

||v + sgn(xi)ei|| 

The additional factor of sgn(xi) in the right-hand side of equation (17.181) assures that 
there is no ambiguity in the sign of the right-hand side of equation (17.151) . In turn, it 
guarantees that all the diagonal elements of the upper-triangular matrix R are positive. 
This is not the definition of Householder refiection used in standard QR decomposition 
routines. Usually, 

i7;(v) = /-2uu*, (7.20) 
with the same u as in (I7.14p . Therefore, 

H'^{^)ir = Tei. (7.21) 

As a consequence, the signs of the diagonal elements of R are random. This is the reason 
why the output of black box QR decomposition routines must be modified in order to 
obtain orthogonal matrices with the correct distribution. 
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Besides being numerically stable, this algorithm has another advantage with respect to 
the Gram-Schmidt orthonormalization. In most applications of numerical analysis Q need 
not be computed explicitly, only Qw does, where w is a specific vector. Generating all the 
Householder refiections is an 0{N'^) process and computing Hj^w requires 0{N) operations 
— it just evaluates the scalar product (u, w). Successively multiplying H^, . . . ,Hi into 
w is an 0{N'^) process. Therefore, it takes in total 0{N'^) operations to compute Qw. 
Instead, the Gram-Schmidt orthonormalization is an 0{N^) process. However, if Q is 
explicitly needed computing the product (17. 7p requires 0{N^) operations too. 

The generalizations to U(A^) and Sp(A^) are straightforward. The only differences are 
in the definitions of the Householder refiections. A suitable choice for U(A^) is 

Hmi^) = -e-'' (/ - 2uu*) . (7.22) 

The unit vector u is 

^ I if) 

V + e ei 



u = II-. II , (7.23) 

||v + e*^ei|| ^ ' 

where v = {xi, . . . , Xm) G C"* and Xi = e*^ The matrix H^iy) is unitary and 



i7„(v)v = ei. (7.24) 

Note that the introduction of e*^ in equations fl7.22p and (17.231) takes into account both the 
potential cancellations and the correct values of the arguments of the diagonal elements 
of the upper-triangular matrix R: equation (I7.24p implies that all the r^jS are real and 
strictly positive. 

For Sp(A^) we have 

H^{w) = -q{I -2uu*). (7.25) 

with 

u = ^, 7.26 

||v + gei|| 

where v = (xi, . . . ^Xm) ^ and Xi = q \\Xi\\. Also in this case 

if^(v)v = ei. (7.27) 



8 A group theoretical interpretation 

We now know how to generate random matrices from any of the classical compact groups 
U(A^), O(A^) and Sp(A^). In order to achieve this goal, we have used little more than 
linear algebra. However simple and convenient this approach is (after all linear algebra 
plays a major role in writing most numerical algorithms), it hides a natural group the- 
oretical structure behind the Householder refiections, which was uncovered by Diaconis 
and Shahshahani [1]. Indeed, generating a random matrix as a product of Householder 
refiections is only one example of a more general method that can be applied to any finite 
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or compact Lie group. Our purpose in this section is to give a flavour of tliis perspective. 
For tlie sake of clarity, as before, we will discuss the orthogonal group 0(iV); the treatment 
of U(A^) and Sp(A^) is, once again, almost identical. 

The need of a more general and elegant approach arises also if one observes that there 
is one feature of the QR decomposition that may not be entirely satisfactory to a pure 
mathematician: Why in order to generate a random point on a N{N — l)/2-dimensional 
manifold — O(A^) in this case — do we need to generate A^^ random numbers? It does 
not look like the most efficient option, even if it is a luxury that can be easily afforded on 
today's computers. 

We will first show how the key ideas that we want to describe apply to finite groups, 
as in this setting they are more transparent. Suppose that we need to generate a random 
element g in a. finite group Tn- In this context, if Fat has p elements, uniform distribution 
simply means that the probability of extracting any g ^ Tn is 1/p. In addition, we assume 
that there exists a chain of subgroups of F^r: 

Fi C F2 C ■■■ C Fat. (8.1) 

In practical situations it is often easier to generate a random element ^ in a smaller sub- 
group, say Tm-i G Fat, than in Ftv itself; we may also know how to take a random repre- 
sentative gm in the left coset Cm = ^m/^m-i- Now, write the decomposition 

Tm — C*m X r^.i. (8-2) 

Once we have chosen a set of representatives of Cm, an element g G F^ is uniquely factorized 
as = gmd, where gm G Cm- If both gm and g are uniformly distributed in Cm and F^-i 
respectively, then g is uniformly distribute in F^. 

We can apply this algorithm iteratively starting from Fi and eventually generate a 
random element in Ftv. In other words, we are given the decomposition 

Fjv = Cjv X ■ • ■ X C2 X Fi. (8.3) 

An element (7 G F^r has a unique representation as a product 

9 = 9n---9u (8.4) 

where gm is a representative in Cm- If the gmS are uniformly distributed in Cm so is g in 
Fat. This is known as the subgroup algorithm [1]. 

This technique applies to random permutations of letters. The chains of subgroups 

is 

{Id} C ^2 C ■ ■ ■ C ^7V, (8.5) 

where Sm is the m-th symmetric group. Other examples include generating random posi- 
tions of Rubik's cube and random elements in GL(A^, Fp), where ¥p is a finite field with p 
elements. 

For O(A^) the decompositions (18.31) and (18. 4p are hidden behind the factorization (17. 7p 
in terms of Householder reflections. Indeed, the subgroup algorithm for O(A^) is contained 
in 
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Theorem 4. Let vi, . . . , v^r be uniformly distributed on Ef^, . . . , E>^ ^ respectively, where 
S"-' = {v„ = (xi, . . . , Xm) e M™| Er=i ^'j = ^} (8.6) 

is the unit sphere in M™. Furthermore, let iJm(v) be the m-th Householder reflection defined 
in equation (I7.18P . The product 

O = H^{y^)Hn^i{vn-i) ■ ■ ■ //2(v2)i/i(vi) (8.7) 

is a random orthogonal matrix with distribution given by Haar measure on 0{N). 

Proof. Suppose we construct O G O(A^) distributed with Haar measure by factorizing a 
matrix X in the Ginibre ensemble as described in section [71 The random matrix O is the 
product of Householder reflections (17.71) and each factor Hmiym) is a function of the unit 
vector Vm G S™'"^ only. We need to show that such v^s are independent and uniformly 
distributed in S"^"^ for m = 1, . . . , A^. 

At each step in the construction of the upper-triangular matrix (17. 6p . the matrix mul- 
tiplied by the m-th Householder reflection, i. e. 

Xrn = Hmi'Vm) " " • H H jy) X , (8.8) 

is still in the Ginibre ensemble. All its elements are, therefore, i.i.d. normal random 
variables. This is a consequence of the invariance 

d^IG{OX) = dficiX), O G 0{N), (8.9) 

of the measure of the Ginibre ensemble. Now, = (xi, . . . ,Xm) is constructed by taking 
the m-th dimensional vector obtained by dropping the first N — m elements of the 
(A^ — m + l)-th column of Xm- The components of are i.i.d. normal random variables. 
It follows that the p.d.f. of is 



1 1 / ™ \ 1 2 



j=l \ j=l J 

Since P{ym) depends only on the length of v^, and not on any angular variable, the unit 
vector = v^/ ||vm|| is uniformly distributed in S*^"^, and is statistically independent 
of Vfc for k ^ m. □ 

Theorem H] is more transparent than relying on the QR decomposition, which seems 
only a clever technical trick. If nothing else, the counting of the number of degrees of 
freedom matches. In fact, the dimension of the unit sphere S'""^ is m — 1. Thus, the total 
number of independent real parameters is 

V^. ^ A^(iV-l) 

J2{m-1)= ^ ^ (8.11) 

m=l 
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Why is theorem m the subgroup algorithm for 0(A^)? As we shall see in theorem [5], 
the factorization (18.71) is unique — provided that we restrict to the definition (17.181) of the 
Householder reflections. This means that 



O(A^) ^§^-1 X ■■■ X §1 X 0(1), 



where 



(8.12) 
(8.13) 

(8.14) 
(8.15) 

(8.16) 



and O G 0(A^ ~ A consequence of theorem H] is that if vjy is uniformly distributed in 
and O is distributed with Haar measure on 0(A^ — 1), then O is Haar distributed 
too. The final link with the subgroup algorithm is given by 

Theorem 5. The left coset space of 0{N — 1) in 0{N) is isomorphic to i.e. 

0{N)/0{N (8.17) 
A complete class of representatives is provided by the map 



0(1) = §° = {-1,1}. 
If we proceed by induction, we obtain 

0{N) = X 0(iV- 1). 

Therefore, a matrix O G O(A^) admits a unique representation as 

O = i/^(v^)n, 

where 

Q = 




i^iv(v) 



— sgn(xi) (/ — 2uu*) if V 7^ ei, 
In if V = ei, 



where 



u 



V + sgn(a;i)ei 
|v + sgn(a;i)ei| 



(8.18) 



(8.19) 



and Xi is the first component of v. 



^The Householder reflections defined in equation (|7.18p are not continuous at ei. Indeed, it can be 
proved that there is no continuous choice of coset representatives. In section [7] this distinction was super- 
fluous: if V is randomly generated, the probability that v — aei is zero. 
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Proof. The group of NxN matrices Q defined in equation f l8.16p is isomorphic to 0(A^ — 1). 
Since 

nei = ei, (8.20) 
0(A^ — 1) can be identified with the subgroup of O(A^) that leave ei invariant, i.e. 

0(A^- 1) = {O G 0(iV)| Oei = ei}. (8.21) 

Now, if two matrices O and O' belong to the same coset, then 

Oei = O'ei = V (8.22) 

and vice versa. In other words, cosets are specified by where ei is mapped. Furthermore, 
since ||Oei|| = 1, we see that they can be identified with the points in the unit sphere. 
Finally, the map fl8.18p is one-to-one and is such that 

HN{v)ei = V. (8.23) 

Therefore, Hj^i spans a complete class of representatives. □ 

Incidentally, theorem H] implies 

Corollary 1. Let (i/io(Af) an-d (i/io(Af-i) be the Haar measures on 0{N) and 0{N — 1) 
respectively. Then 

dfJ-oiN) = dfigN-i X (i/io(Ar-i), (8.24) 
where dfigN-i is the uniform measure on E>^^^. 

What is the meaning of dfi§N-i7 Given that we are dealing with uniform random 
variables, it is quite natural that we end up with the uniform measure. In this case, 
however, it has a precise group theoretical interpretation. Left multiplication of the right- 
hand side of equation fl8.15p by O' G O(A^) induces a map on the coset space: 

0'HN{v)n = HN{^')n'n = HN{v)n". (8.25) 

Since the decomposition fl8.15p is unique the transformation v i-^ v' is well defined. This 
map can be easily determined. A coset is specified by where ei is mapped, therefore 

0'ifjv(v)ei = O'v = V = i/Ar(v')ei. (8.26) 

If V is uniformly distributed on the unit circle so is v' = Ov. Thus, dfi^N-i is the unique 
measure on the coset space 0{N)/0{N — 1) invariant under the left action of O(A^). Its 
uniqueness follows from that of Haar measure and from the factorization 08.240 . 

Corollary [T] is a particular case of a theorem that holds under general hypotheses for 
topological compact groups. Indeed, let F be such a group, S a closed subgroup and 
C = F/S. Furthermore, let d/xr, dfj,c and dfix be the respective invariant measures, then 

dfir = djji-E X due- (8.27) 
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